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Abstract 


A  novel  correlative  technique  is  introduced  for  comparison  of  measurements  of  similar 
quantities  made  using  different  techniques.  This  process  involves  a  generalized 
least-squares  fitting  method  which  can  be  used  to  estimate  the  slope  of  the  best-fitting 
straight  line  that  results  when  two  separate  data  sets  that  are  expected  to  be  linearly 
correlated  are  compared  via  scatter-plots.  The  different  data  types  may  be  subject  to 
different  uncertainties  in  their  measurements.  The  technique  determines  and  graphs  the 
relationship  between  the  errors  in  each  method  and  the  slope  of  the  line  of  best  fit, 
assuming  Gaussian  statistics.  We  hope  that  this  new  technique  will  find  appUcation  in  a 
wide  range  of  situations,  e.g.,  ranging  from  radar  and  satelUte  measurements  to  any 
physical  measurements.  The  technique  can  also  be  used  to  calibrate  multi -sensor,  radar 
and  satelUte,  systems. 
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Resume 


Une  technique  correlative  innovatrice  de  comparaison  de  mesures  de  quantites 
similaires  par  des  methodes  differentes  est  presentee.  Ce  processus  fait  intervenir  une 
methode  generalisee  d’ajustement  par  les  moindres  carres  pour  Testimation  de  la  pente 
de  la  droite  la  mieux  ajustfe  resultant  de  la  comparaison  par  diagrammes  de  dispersion 
de  deux  jeux  de  donnees  distincts  dont  on  s’attend  a  ce  qu’ils  presentent  une  correlation 
lineaire.  Les  mesures  permettant  d’obtenir  differents  types  de  donnees  peuvent  etre 
sujettes  a  differentes  incertitudes.  La  technique  proposee  permet  de  determiner  et  de 
representer  graphiquement  la  relation  entre  les  erreurs  pour  chaque  methode  et  la  pente 
de  la  droite  la  mieux  ajustee  en  supposant  une  repartition  stochastique  gaussienne. 

Nous  esperons  que  cette  nouvelle  technique  sera  applicable  a  une  gamme  etendue  de 
situations,  p.  ex.  depuis  les  mesures  effectuees  au  moyen  du  radar  et  de  satellites 
jusqu’a  toute  mesure  physique.  Cette  technique  peut  egalement  servir  a  Tetalonnage  de 
systemes  multicapteurs,  radar  et  satellites. 


ii 
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Executive  summary 


One  of  the  problems  that  frequently  confronts  an  experimentalist  is  that  of  fitting  a 
straight  line  to  data,  that  is,  the  problem  of  determining  the  functional  relationship 
between  two  variables.  Least-squares  techniques  for  examining  best-fit  lines  to 
correlated  data  sets  are  well  established,  but  in  general  are  restricted  to  cases  in  which 
the  errors  in  the  two  data  sets  being  compared  are  known.  The  simplest  case  is  that  in 
which  there  is  no  error  in  one  variable  (usually  plotted  as  the  abscissa),  and  the  error  in 
the  other  variable  is  known.  This  procedure  is  described  in  almost  all  elementary  books 
on  statistics. 

A  more  complex  case  is  that  in  which  both  variables  have  errors,  but  both  errors  are 
well  known.  Examples  which  demonstrate  how  to  deal  with  such  cases  have  been 
shown  in  the  literature.  However,  there  are  times  when  two  data  sets  are  compared  in 
which  the  intrinsic  measurement  error  of  one  or  both  variables  are  unknown.  In  this 
case,  there  is  no  algorithm  available  to  determine  the  best-fit  line.  The  situation  can  be 
further  clouded  if  the  different  techniques  measure  similar,  but  not  identical, 
parameters.  Our  purpose  in  this  report  is  to  consider  optimal  ways  to  compare  such 
data. 

A  novel  correlative  technique  is  introduced  for  comparison  of  measurements  of  similar 
quantities  made  using  different  techniques.  This  process  involves  a  generaUzed 
least-squares  fitting  method  which  can  be  used  to  estimate  the  slope  of  the  best-fitting 
straight  line  that  results  when  two  separate  data  sets  which  are  expected  to  be  linearly 
correlated  are  compared  via  scatter-plots.  The  different  data  types  may  be  subject  to 
different  uncertainties  in  (heir  measurements.  The  technique  determines  and  graphs  the 
relationship  between  the  errors  in  each  method  and  the  slope  of  the  line  of  best  fit, 
assuming  Gaussian  statistics.  We  hope  that  this  new  method  will  find  application  in  a 
wide  range  of  situations,  e.g.,  ranging  from  radar  and  satellite  measurements  to  any 
physical  measurements.  The  technique  can  also  be  used  to  calibrate  multi-sensor,  radar 
and  satellite,  systems. 


T.  Thayaparan,  W.  K.  Hocking,  S.  J.  Franke.  2001.  A  Novel  Method  for  Statistical 
Comparison  of  Geophysical  Data  by  Multiple  Instruments  Which  Have  Differing  Accuracies. 
DREO  TM  2001-147.  Defence  Research  Establishment  Ottawa. 
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Sommaire 


L’ajustement  d’une  droite  a  des  donnees,  c’est’a-dire  la  determination  de  la  relation 
fonctionnelle  entre  deux  variables,  constitue  Tun  des  problemes  auxquels  sont 
frequemment  confrontes  les  experimentateurs.  Les  methodes  d’examen  par  les 
moindres  carres  des  droites  les  mieux  ajustees  a  des  jeux  de  donnees  sont  couramment 
utilisees,  mais  en  general  d’application  limitee  aux  cas  dans  lesquels  les  erreurs  dans 
les  deux  jeux  de  donnees  comparees  sont  connues.  Le  cas  le  plus  simple  est  celui  dans 
lequel  une  des  variables  ne  comporte  pas  d’erreur  (habituellement  celle  qui  est 
representee  en  abscisse)  et  Terreur  entachant  Tautre  variable  est  connue.  Cette 
procedure  est  decrite  dans  presque  tous  les  ouvrages  elementaires  de  statistique. 

La  situation  est  plus  complexe  lorsque  des  erreurs  entachent  les  deux  variables,  mais 
que  ces  erreurs  sont  bien  connues.  Des  exemples  de  la  maniere  dont  il  convient  de 
traiter  ces  situations  ont  ete  documentes.  Cependant,  dans  certains  cas,  il  faut  comparer 
deux  jeux  de  donnees  pour  lesquels  I’erreur  intrinseque  de  mesure  d’une  ou  des  deux 
variables  n’est  pas  connue.  Dans  cette  situation,  il  n’existe  aucun  algorithme  permettant 
de  determiner  la  droite  de  meilleur  ajustement.  Le  probleme  devient  encore  plus 
epineux  si  des  methodes  dilferentes  sont  appliquees  pour  mesurer  des  parametres 
similaires  mais  non  identiques.  Le  present  rapport  a  comme  objet  Texamen  de 
methodes  optimales  de  comparaison  de  telles  donnees. 

Une  technique  correlative  innovatrice  de  comparaison  de  mesures  de  quantites 
similaires  par  des  methodes  differentes  est  presentee.  Ce  processus  fait  intervenir  une 
methode  generalisee  d’ajustement  par  les  moindres  carres  pour  I’estimation  de  la  pente 
de  la  droite  la  mieux  ajustee  resultant  de  la  comparaison  par  diagrammes  de  dispersion 
de  deux  jeux  de  donnees  distincts  dont  on  s’attend  a  ce  qu’ils  presentent  une  correlation 
lineaire.  Les  mesures  permettant  d’obtenir  differents  types  de  donnees  peuvent  etre 
sujettes  a  differentes  incertitudes.  La  technique  proposee  permet  de  determiner  et  de 
representer  graphiquement  la  relation  entre  les  erreurs  pour  chaque  methode  et  la  pente 
de  la  droite  la  mieux  ajustee  en  supposant  une  repartition  stochastique  gaussienne. 

Nous  esperons  que  cette  nouvelle  technique  sera  applicable  a  une  gamme  etendue  de 
situations,  p.  ex.  depuis  les  mesures  effectuees  au  moyen  du  radar  et  de  satellites 
jusqu’a  toute  mesure  physique.  Cette  technique  peut  egalement  servir  a  I’etalonnage  de 
systemes  multicapteurs,  radar  et  satellites. 


T.  Thayaparan,  W.  K.  Hocking,  S.  J.  Franke.  2001.  Une  methode  innovatrice  de  comparaison 
statistique  de  jeux  de  donnees  geophysiques  recueillis  au  moyen  d’instruments  offrant 
differentes  exactitudes.  DREO  TM  2001--147.  Centre  de  Recherches  pour  la  Defense  Ottawa. 
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1 .  Introduction 


Least-squares  techniques  for  examining  best-fit  lines  to  correlated  data  sets  are  well 
established,  but  in  general  are  restricted  to  cases  in  which  the  errors  in  the  two  data  sets 
being  compared  are  known.  The  simplest  case  is  that  in  which  there  is  no  error  in  one 
variable  (usually  plotted  as  the  abscissa),  and  the  error  in  the  other  variable  is  known. 
This  procedure  is  described  in  almost  all  elementary  books  on  statistics  [1-5]. 

A  more  complex  case  is  that  in  which  both  variables  have  errors,  but  both  errors  are 
well  known.  Examples  which  demonstrate  how  to  deal  with  such  cases  have  been 
shown  by  (amongst  others)  by  [6-16]. 

However,  there  are  times  when  two  data  sets  are  compared  in  which  the  intrinsic 
measurement  error  of  one  or  both  variables  are  unknown.  In  this  case,  there  is  no 
algorithm  available  to  determine  the  best-fit  line.  The  situation  can  be  further  clouded  if 
the  different  techniques  measure  similar,  but  not  identical,  parameters.  An  example  is 
shown  in  Figure  1,  where  two  time-series  of  the  eastward  wind  component  in  the 
atmosphere  at  an  altitude  of  85  km  are  shown.  These  data  were  recorded  using  two 
different  techniques,  one  being  a  spaced  antenna  procedure  applied  to  radiowave 
diffraction  patterns  produced  by  MF  (medium  frequency)  radio  scatter  from 
ionospheric  irregularities  [17],  and  the  other  using  interferometric  applications  applied 
to  meteor  trails  formed  in  the  Earth’s  atmosphere  (e.g.,  see  [18]  and  references 
there-in).  Although  each  procedure  ostensibly  measures  “winds"  at  85  km,  the  two 
instruments  measure  over  different  areas  and  depths  of  the  sky.  As  a  result,  there  is  a 
possibility  that  they  measure  different  phenomena,  especially  if  the  winds  show 
significant  spatial  variation.  Our  purpose  in  this  paper  is  to  consider  optimal  ways  to 
compare  such  data.  It  should  be  emphasized  here  that  the  new  proposed  technique  is 
not  restricted  to  above  example  but  can  be  applied  in  a  wide  range  of  situations,  e.g., 
ranging  from  radar  and  satellite  measurements  to  any  physical  measurements.  The 
technique  can  also  be  used  to  calibrate  multi-sensor,  radar  and  satellite,  systems. 

This  report  begins  by  developing  the  required  nomenclature,  and  a  theoretical  basis  for 
the  calculations.  The  basic  theory  for  the  ensuing  procedures  are  then  developed.  In 
order  to  confirm  the  proposed  theory,  we  simulate  data  using  computer  techniques,  in 
which  the  intrinsic  errors  are  pre-specified.  We  demonstrate  the  effects  of  employing  at 
least  one  existing  method  of  multi-error  regression  analysis  [6]  to  these  data  sets,  and 
point  out  the  limitations  of  such  procedures.  Finally,  we  apply  the  proposed  procedures 
to  two  sample  data  sets,  including  the  points  shown  in  Figure  1. 
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Zonal  Wind  Time  Series  Comparison,  85  km  altitude. 

MF  winds  vs  Meteor  Winds. 

London,  Ontario,  Canada,  June  10-18,  1999. 


Day  in  June  1999 


Figure  1:  An  example  of  two  time  series  of  data  obtained  by  two  different  models. 
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2.  Nomenclature  and  Theory 


The  goal  in  this  section  is  to  compare  two  data  sets  that  are  assumed  to  be  linearly 
related,  and  to  make  meaningful  statements  about  the  relationship  between  them 
without  having  any  a-priori  knowledge  about  system  errors.  For  the  purposes  of  our 
theoretical  developments,  we  will  denote  these  data  sets  as  {a:*}  and  {yj}  (*  =  L.iV). 
The  quantities  {xj}  are  assumed  to  represent  measurements  of  a  randomly  varying 
parameter,  denoted  here  by  {vj},  and  the  quantities  {yj}  represent  measurements  of  a 
parameter  which  is  assumed  to  be  linearly  related  to  The  data  points  {y,}  are  in 
fact  assumed  to  have  a  constant,  but  unknown,  gain  relative  to  the  {vi},  which  we  will 
denoted  as  yo-  We  will  assume  that  each  data  set  has  a  mean  of  zero:  in  practice,  if  this 
is  not  true,  it  is  easily  rendered  so  by  removing  the  means  from  the  data  points.  The 
model  equations  are  therefore 


(1) 


Xi  =  Vi  +  Sxi 


(2) 


pi  =  go  Vi  +  6yi 


where  6xi  and  &yi  are  random  (noise)  components.  We  shall  use  the  notation  W(y,  a^) 
to  denote  the  Gaussian  probability  density  function  with  mean  y  and  variance  cr^.  The 
following  equations  summarize  our  assumptions  about  the  distribution  of  the 
measurement  errors  and  the  intrinsic  associated  parameters: 


(3) 


Vi  =  N{Q,^l) 


(4) 


6xi  = 


(5) 


%  =iV(0,cr^) 


We  also  assume  that  Vi,  Sxi  and  Spi  are  mutually  independent.  Our  purpose  now  is  to 
determine  an  expression  relating  the  variable  go,  (tI  and 

We  now  square  Eq.  1  and  Eq.  2,  and  take  ensemble  averages  (<>)  to  produce  the 
following  equations  (remembering  that  we  have  assumed  zero-mean  quantities): 
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(6)  <  a:  ■  >  =  <  v?  >  +  <  6Xi  >  = 

(7)  <y^>=  gl<v^>  +  <6y1>=  glYil  +  al 


Furthermore,  if  we  also  multiply  Eq.  1  by  Eq.  2,  and  then  take  ensemble  averages,  we 
produced  the  additional  equations 


(8) 


<  Xipi  >=  go  <  Vi  >=  go 


Note  that  we  have  assumed  that  Vi,  6xi  and  6yi  are  uncorrelated,  so  the  averages  of 
crossed  terms  are  zero. 

In  what  follows,  the  ensemble-averaged  quantities  involving  the  measured  data  are 
replaced  with  sample  expectations.  For  example,  <  x?  >  is  replaced  with 

=  {EiXi)/N,  <  yl  >  is  replaced  with  =  (Si2/?)/JV  and  <  Xiyi  >  is  replaced 
with  =  {T,iXiyi)IN. 

We  also  introduced  two  more  variables  in  the  following  way.  We  recognize  that  the 
quantities  are,  in  essence,  the  values  of  {xi}  prior  to  the  addition  of  the  random 
(noise)  component  and  likewise  the  values  {goVi}  are  the  values  of  {^i}  prior  to 

the  addition  of  the  noise  component  We  may  consider  these  “non-noise" 
quantities  as  the  “signal"  embedded  in  our  data.  We  henceforth  use  the  definitions 

(9)  El  = 


(10) 


where  and  represent  the  variances  of  the  “signal"  component  of  the  {^i}  and 
{j/i}  respectively.  Then  Eq.  6  and  Eq.  7,  combined  with  our  definitions  of  ^  and 
produce  the  result 


(11) 


2 

X 
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and 


(12)  ey  =  ^  +  <4 


i.e.,  they  tell  us  the  somewhat  obvious  result  that  the  overall  variance  of  each  data  set  is 
the  sum  of  the  variance  associated  with  the  signal  and  variance  associated  with  the 
noise.  Likewise  we  also  have  the  result 


(13) 


^xy  9o^x 


We  now  turn  to  regression  estimators.  In  many  studies  of  regression,  it  is  common  to 
treat  one  variable  as  if  it  has  no  error,  and  assume  all  the  error  is  associated  with  the 
second  variable.  For  example,  one  might  assume  that  era?  =  0  in  our  above  equations, 
while  letting  any  error  be  associated  with  the  variable  yi.  Whilst  this  is  untrue  of  our 
data,  it  makes  a  convenient  starting  point.  If  we  take  Eq.  1 1  and  Eq.  13  and  assume  that 
ax  —  0,  and  solve  for  go,  then  we  produce  an  estimator  for  the  slope  which  we  will  call 
Qx  .  This  quantity  is  given  by 


(14) 


9x  —  Cxy/ C 


2 

X 


This  is  exactly  the  same  parameter  as  that  deduced  in  any  standard  text  book  as  the 
“slope  of  the  leasLsquares  best  fit  line"  i.e.,  the  regression  of  y  on  x  [2,5].  The  “best-fit" 
line  has  an  equation 


(15) 


9x  “b  ^1 7 


where  a  is  a  constant. 

Likewise,  by  assuming  that  ay  =  0,  and  solving  for  go  in  Eq.  12  and  Eq.  13,  we 
produce  the  slope  corresponding  to  the  regression  of  x  on  y,  which  obeys  the  equation 


(16) 


=  {i/9y')yi  +  C2, 


viz. 
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Note  that  in  the  second  case  we  have  used  an  inverted  form  for  the  slope  (i.e.,  it 
involves  l/gy),  so  that  we  can  write 


(18) 


Vi  —  Py  ^3? 


Whilst  neither  of  these  estimators  represent  the  true  slope  for  our  situation,  they  make 
excellent  starting  points  for  a  derivation  of  a  better  estimator  for  g^.  In  all  real  cases, 
gj  is  less  than  or  equal  to  and  gy  is  greater  than  or  equal  to  po-  The  values  g^  and 
gy  are  also  very  easy  to  find,  since  many  software  packages  produce  these  quantities  as 
standard  output  (recognizing  that  the  regression  of  x  on  y  needs  to  have  its  slope 
inverted  to  be  consistent  with  our  definition  of  gy).  We  emphasize  that,  although  these 
quantities  were  derived  assuming  either  that  ajc  =  0  or  that  Uy  =  0,  we  will  not  assume 
this  in  our  more  general  analysis.  These  two  slope  estimators  represent  only  a  starting 
point  for  our  calculations.  For  any  data  set,  gx  and  gy  are  only  estimators  of 
population  values  gx  {=<  gj  >)  and  gy  {—<  gy  >),  but  as  long  as  N  is  moderately 
large,  we  can  take  our  estimates  to  be  reasonable  approximations  of  the  population 
values,  in  the  same  manner  that  is  usually  applied  to  statistical  studies. 

Eq.  14  tells  us  that  ^  Eq.  3  gives  =  go^l^  so  we  may  equate  these 

two  equations  to  give  =  poS|.  Applying  Eq.  1 1  and  eliminating  then  gives 


(19) 


^x/^x  —  (po/px  1)^^  • 


Likewise, 


(20) 


^vl^y  —  (py/po  !■)  ^  • 


Thus  we  see  that  the  ratios  of  the  “noise"  variance  to  the  “signal"  variance  are  simple 
functions  of  50  gx  for  iho  x-series,  and  simple  functions  of  go  and  gy  for  the 
y- series.  Furthermore,  we  can  use  Eq.  1 1  and  Eq.  12  to  write 

(21)  a^  =  ex(l-5x/po)'/'- 


DREOTM  2001-1 47 


Similarly, 


(22) 


<^1/ =  ^1/(1  - 


Since  and  are  just  sample  variances,  and  estimators  for  and  gy  can  easily  be 
determined  for  any  data-set,  the  only  unknown  quantities  in  the  above  two  equations  are 
go,  (Tx  and  ay.  Since  these  three  quantities  are  uniquely  inter-related;  specification  of 
any  one  of  them  allows  the  determination  of  the  other  two.  This  is  the  key  point  of  our 
procedure,  and  is  a  feature  which  we  will  employ  throughout  our  future  experimental 
examples. 
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3.  Monte  Carlo  Simulations 


In  order  to  confirm  our  calculations  above,  we  now  turn  to  Monte  Carlo  modelling.  Our 
simulations  begin  by  developing  a  data-set  of  points  {vi  }  which  have  a  Gaussian 
probability  distribution  with  zero  mean,  as  assumed  in  our  previous  theory  We  assume 
that  a  Gaussian  probability  distribution  is  suitable  for  our  data.  The  standard  deviation 
of  the  set  is  written  as  Then,  for  each  point  Vi  ,  we  calculate  a  point  Vgi  —  QqVi  ,  so 
that  the  set  {vgi}  has  zero  mean,  and  represents  a  re-scaling  of  the  original  set  by 
time.  Consistent  with  our  description  in  the  previous  section,  we  consider  these  sets  of 
points  to  be  {xi}  and  {y/}  values  prior  to  the  addition  of  random  fluctuations,  and 
therefore  denote  the  standard  deviation  of  the  original  set  {vi  }  as  and  the  standard 
deviation  for  the  quantities  {vg/}  as  =  go^v  =  These  standard  deviations 
essentially  represent  the  “natural”  spread  of  the  (simulated)  data  set  which  we  are 
generating.  At  this  stage  all  the  points  ,  yi  }  lie  in  a  straight  line  if  plotted  in  a 
“scatter  plot”  format.  We  now  add  a  random  component,  by  adding  a  vector  Syi) 
to  each  point,  where  Sxi  is  a  randomly  selected  value  from  a  data-set  with  a  Gaussian 
probability  distribution  of  zero  mean  and  standard  deviation  and  similarly  Syi  is  a 
randomly  selected  value  from  a  separate  data-set  with  a  Gaussian  probability 
distribution  of  zero  mean  and  standard  deviation  ay.  This  produces  scatter  plots  which 
look  like  those  shown  in  Figure  2.  We  will  denote  this  new  set  as  {xi,yi},  consistent 
with  the  above  theory. 

Following  this,  we  perform  two  regression  analyses.  First  we  treat  the  {xi}  as  if  they 
have  no  errors,  and  produce  a  best-fit  line  with  slope  gx  (viz.  the  equation  of  this  line  is 
equation  yi  —  gxXi  +  Ci).  This  is  the  result  of  the  regression  of  y  on  x.  We  then  treat 
the  {yi}  points  as  if  they  have  no  error,  and  produce  a  best-fit  line  of  slope  gy  (i.e.,  we 
find  the  best-fit  line  satisfying  xi  —  g^yi  -h  C2,  but  then  convert  it  to  the  form 
yi  =  gyXi  +  C3,  where  gy  =  1/52)-  Examples  are  shown  in  Figure  2.  It  should 
especially  be  noted  that  the  regression  of  y  on  a:  produces  a  slope  which  is  less  than  the 
true  slope,  and  likewise  the  regression  of  x  on  y  also  produces  an  underestimate  of  the 
true  slope  i.e.,  and  g^  <  yo-  In  the  second  case  y2  <  yo  means  that  y^  >  yo-  In 

addition  to  performing  these  simulations,  we  also  applied  some  existing  dual-error 
algorithms  to  our  data.  In  particular,  we  have  applied  the  procedure  of  [6],  using  the 
known  pre-assigned  errors  for  each  data-set. 

We  have  performed  the  above  calculation  on  many  thousands  of  sets  of  points,  and  have 
determined  that  values  of  y^  and  gy  for  different  combinations  of  Sy,  nnd  ay. 
Figure  3  summarizes  our  results.  It  shows  not  only  the  average  values  of  gx  and  gy  for 
many  different  realizations  for  the  case  yo  —  1  and  pre-specified  values  of  ax  and  ay 
(Figures  3a  and  3c),  but  also  the  standard  deviations  (i.e.,  spread)  in  each  estimate  of  y^? 
and  gy  (right  hand  graphs).  The  results  of  the  fitting  procedures  of  [6]  are  also  shown  as 
Figures  3e  and  3f.  Figures  3d  and  3f  serve  to  demonstrate  that  [6]’s  procedure  works 
well,  in  that  the  slopes  are  always  close  to  unity,  but  the  method  did  require  prior 
knowledge  about  the  errors  associated  with  each  of  the  abscissa  and  the  ordinate.  If 
these  errors  are  not  known,  the  method  cannot  be  applied,  and  our  intent  here  is  to  deal 
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Figure  2:  Sample  simulated  scatter  plots.  In  each  case,  we  have  taken  E®  =  Ey  =  40,  and  used  values  forax  anof 
(Ty  as  shown.  The  slopes  Qx  and  py ,  which  were  determined  by  standard  least-squares  regression  analysis,  are 
indicated  on  the  figure.  The  slope  gxy  is  the  result  of  a  fit  due  to  York  [6],  which  required  prior  knowledge  of  the 
values  forcTx  and<Ty. 
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Figure  3:  Mean  values  of  (a)  gx  and  (c)  gy  calculated  from  our  simulations,  as  functions  of  a  x  and  <7y,  where  in 
each  case  the  means  were  determined  from  several  hundred  realizations  (each  realization  being  like  the  example 
shown  in  Figure  2,  but  each  time  using  different  sets  of  random  numbers).  We  have  usedllx  =  Ey  =  40.  In  all 
cases  the  abscissa  is  rjx  and  the  ordinate  is  (7y,  and  the  values  range  from  0  to  80.  Figures  (b)  and  (d)  show  the 
standard  deviations  determined  during  these  simulations,  where  (b)  shows  the  standard  deviations  corresponding  to 
(a),  and  (d)  shows  the  standard  deviations  for  (c).  Figure  (e)  shows  the  results  of  calculations  of  the  slope  of  the 
best-fit  line  using  the  procedure  of  York  [6],  where  the  known  errors  have  been  used  in  the  calculations,  and  (f) 
shows  the  associated  standard  deviations  for  this  procedure. 
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Figure  4:  Comparison  of  our  calculated  values  of  a  xj'^x  as  a  function  of  Qx  I  go,  compared  to  our  theoretical 
determinations  (Eq  19). 


with  situations  with  unknown  errors.  We  shall  therefore  pursue  the  procedures  of  [6]  no 
further,  and  will  turn  to  a  deeper  discussion  of  Figures  3a,  3b,  3c  and  3d. 

With  regard  to  Figures  3a  and  3c,  it  should  be  noticed  in  particular  that  the  contours  are 
vertical  in  the  first  case,  indicating  that  gx  is  a  function  only  of  for  a  specified  value 
of  go,  Sjc  and  Ey,  (which  is  consistent  with  Eq.  19),  and  the  contours  are  horizontal  for 
the  second  case,  indicating  that  gy  is  independent  of  ay  (consistent  with  Eq.  20).  We 
have  also  been  able  to  show  that  a  graph  of  a^j'^x  as  a  function  of  par/Po  fias  a  form 
like  that  shown  in  Figure  4,  irrespective  of  go.  Ex  and  ax-  Similarly  we  may  determine 
graphs  of  ayjlly  as  a  function  of  po/py  As  shown  in  Figure  4,  our  simulations  match 
very  closely  Eq.  19  and  Eq.  20.  Figure  4,  and  Eq.  19  and  Eq.  20  are  the  key  to  our 
calculations  performed  in  this  paper.  We  now  turn  to  the  application  of  these  formulae 
in  real  situations. 

To  illustrated  our  applications  to  experimental  measurements,  we  turn  to  Figure  5, 
which  shows  a  scatter  plot  for  an  extended  data  series  which  includes  the  points  in 
Figure  1 .  The  data  include  both  meridional  and  zonal  components  of  the  wind 
measured  by  co-located  meteor  and  MF  (Medium  Frequency)  radars  at  London, 

^  Ontario,  Canada  (see  [18]  for  specific  details  about  these  innstruments).  The  graph  in 

Figure  1  shows  only  the  zonal  component  for  a  short  period  of  8  days,  but  our 
comparison  uses  a  super-set  of  these  data,  covering  the  entire  period  from  June  1  to 
July  31,  1999  (with  some  data  gaps  during  which  the  radars  were  used  for  other 
purposes).  For  the  present  we  will  not  concern  ourselves  with  the  physical  meaning  of 
the  data,  but  will  simply  regard  them  as  representative  of  any  such  correlative  study 
between  two  such  data  sets.  However,  we  do  point  out  that  the  winds  at  85  km  altitude 
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Figure  5:  Results  of  applying  the  fitting  procedure  to  the  data  like  those  shown  in  fig.  1.  The  left-hand  graph  (a) 
shows  the  scatter  plot  of  MF  wind  components  vs.  meteor  wind  components  for  an  altitude  of  85  km  for  the  period  of 
June  and  July,  1999,  using  co-hcated  radars  near  London  Ontario,  Canada.  Relevant  parameters  (offsets,  g-x  and 
gy)  are  also  shown  on  the  figure.  The  quantities  xa  andya  refer  to  axis  intercepts,  as  described  in  the  text.  Errors 
are  quoted  at  the  two-sigma  level  (approximately  95%  significance).  Fig  (b)  shows  the  relationship  between  go,  <7 x, 
and  <7 y,  as  described  in  the  text. 


in  the  upper  atmosphere  are  highly  variable,  both  temporally  and  spatially,  and 
comprise  multiple  large-amplitude  waves.  Vertical  excursions  as  large  as  5-6  km  in  a 
few  hours  are  not  uncommon.  Hence  the  reader  should  not  be  surprised  by  an  apparent 
low  correlation  between  our  two  data  sets;  the  higher  correlations  which  are  typical  of 
(for  example)  tropospheric  data  are  simply  not  possible  in  upper  atmospheric  cases. 

When  presented  with  real  data,  we  perform  the  following  calculations.  We  first  find  the 
standard  deviation  of  each  dataset  {iCi}  and  {j/i},  which  we  will  denote  as  and  ^y. 
We  then  also  determine  the  slopes  of  the  regression  of  y  on  x,  to  give  Qx,  and  the 
regression  of  x  on  y,  thereby  obtaining  Qy  (remembering  that  the  slope  of  the  regression 
of  a;  on  2/  must  be  inverted  to  give  gy).  We  also  determine  intercepts  of  these  best-fit 
lines  with  the  axes,  which  are  denoted  as  Xoi  and  t/oi  in  this  (and  subsequent)  figures. 
The  value  xoi  refers  to  the  intercept  of  the  line  of  slope  Qx  with  the  ordinate,  and  the 
value  2/oi  refers  to  the  intercept  of  the  line  of  slope  Qy  with  the  abscissa.  Following  this, 
we  assume  various  values  for  ctx,  starting  at  0.0  and  stepping  up  in  small  steps.  For 
each  value  of  cr^,  we  calculate  a  value  for  Ex  through  the  relation  Ex  —  {^x  ~ 

(from  Eq.  11).  We  can  then  deduce  the  ratio  (Txf'^x^  and  from  the  graph  in  Figure  3,  or 
from  Eq.  19,  we  therefore  deduce  Qx/go-  Knowing  gx,  we  can  therefore  deduce  go  for 
this  particular  choice  of  Following  this,  we  can  deduce  Ey  because  Ey  ~ 

Finally,  we  may  deduce  ay  —  {^y  —  We  may  therefore  plot  the  relationship 

between  go,  ax  and  ay  as  shown  in  Figure  5b.  Note  that  the  upper  axis  {ay)  requires  a 
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Figure  6:  As  for  Figure  5,  but  using  data  from  a  tropospheric  wind-profiler  (see  text). 


different  scaling  for  each  new  data-set  for  which  this  procedure  is  applied. 

This  procedure  therefore  explains  how  the  secondary  graphs  in  Figure  5  were  produced. 
Extra  information  is  required  in  order  to  determine  the  “correct"  combination  of  go,  ax 
and  ay.  (Thayaparan  and  Hocking  [19]  attempted  to  use  an  additional  graph  like  that 
shown  in  Figure  4,  but  relating  ayjlhy  to  to  try  and  uniquely  determine  go,  ax 
and  ay,  and  initial  trials  suggested  some  success.  However,  it  subsequently  evolved  that 
the  information  available  with  this  process  was  redundant,  and  a  unique  determination 
was  impossible  without  further  information.  Therefore  the  reader  should  be  aware  that 
Figure  4  of  [19],  is  incorrect,  as  are  the  developments  following  Eq.  9  of  that  work.) 

In  all  of  our  scatter  plots  we  will  henceforth  include  graphs  like  that  shown  in  Figure 
5b.  This  graph  encompass  many  of  the  features  of  other  fitting  algorithms.  For 
example,  if  the  user  has  sufficient  justification  to  assume  ax  =  ay,  then  can  be  read 
from  the  graph  directly.  If  a  particular  value  of  ax  (or  for  that  matter  ay)  is  known,  then 
again  go  can  be  found  uniquely.  If  it  is  known  that  go  is  1.0  (for  example),  then  ax  and 
ay  can  be  immediately  determined.  The  so-called  “OLE”  (Ordinary  Least-squares 
Bisector”)  method  (in  which  the  sum  of  the  squares  of  the  perpendicular  distances 
between  the  best-fit  line  and  the  points  is  minimized)  is  another  subset  of  this  curve, 
and  corresponds  to  the  case  where  the  ratio  of  ax  to  ay  is  assumed  to  be  proportional  to 
po-  In  Figure  5b,  the  cases  where  ax  =  ay  and  go  =  1^0  are  almost  coincident,  but  this 
need  not  always  be  so.  Thus  this  graph  encompasses  many  of  the  features  available 
with  other  least  squares  algorithms,  but  gives  a  much  clearer  picture  of  the  relationship 
between  these  various  variables.  Most  other  algorithms  discussed  in  the  current 
literature  produce  subsets  of  the  more  general  {go,  ax,  ay)  combinations  shown  in 
figures  like  Figure  5b. 
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Figure  6  shows  another  example  of  application  of  this  process.  In  this  case,  the  data 
were  taken  using  a  tropospheric  radar,  so  the  scatter  is  less  severe.  This  is  true  because 
the  wind  fields  in  the  troposphere  tend  to  be  much  less  chaotic  than  those  in  the  upper 
atmosphere.  The  data  refer  to  radial  Doppler  velocities  measured  on  two  beams 
pointing  at  10.9  degrees  off-vertical  in  opposite  directions.  Again,  the  various  possible 
combinations  of  and  Cy  are  defined  by  Figure  6b.  In  this  case,  because  the  beams 
are  essentially  identical,  it  would  be  sensible  to  surmise  that  ax  —  ay,  so  that  we  may 
determine  that  po  =  —1.028.  This  therefore  tells  us  that  the  effective  beams  in  each 
case  are  not  at  the  same  zenithal  angle,  but  differ  by  about  18  minutes  of  arc.  The  exact 
reason  for  this  is  not  important  for  this  paper,  but  our  statistical  procedure  has  been  able 
to  determine  this  asymmetry. 
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4.  The  Meaning  of  Ux  and  ay _ 

Finally,  we  need  to  consider  what  ax  and  ay  really  represent.  If  we  are  using  two 
instruments  to  measure  the  same  quantity,  then  these  two  variables  tell  us  the 
measurement  error  in  each  technique.  However,  in  many  situations  in  geophysics  we 
5  cannot  measure  exactly  the  same  phenomenon.  The  comparison  between  the  MF  and 

meteor  winds  which  we  considered  in  Figure  5  is  one  such  example.  Whilst  it  would  be 
desirable  for  each  instmment  to  measure  in  exactly  the  same  region  of  space,  over  the 
^  same  interval  of  time,  this  is  simply  not  possible.  Thus  the  values  of  the  quantities  ax 

and  ay  reflect  not  only  the  intrinsic  measurement  errors  of  each  instrument,  but  also 
contain  information  about  how  “different"  the  two  techniques  are.  In  essence,  the 
least-squares  fitting  procedure  indirectly  assumes  that  there  is  a  “compromise" 
measurement  technique  that  is  neither  the  MF  radar  nor  the  meteor  method,  and 
assumes  that  each  instrument  is  trying  to  achieve  this  compromise.  For  example,  as  a 
crude  example,  the  MF  radar  measures  in  a  volume  of  space  which  is  typically  10-20 
km  across,  whilst  the  meteor  radar  measures  over  a  region  of  space  of  perhaps  a  100  km 
or  more  in  diameter.  This  “compromise"  measurement  might  represent  an  instrument 
which  could  measure  in  an  intermediate  volume  of  space  -  perhaps  with  a  width  of  fifty 
km  or  so  and  a  depth  of  3-4  km.  Of  course  this  grossly  oversimplifies  the  true  situation, 
but  illustrates  the  concept.  In  reality,  we  can  never  actually  state  what  this  “compromise 
technique"  is,  or  even  achieve  it.  However,  we  can  say  that  ax  and  ay  contain 
information  about  both  the  intrinsic  measurement  errors  of  each  technique  and  also  the 
level  of  departure  of  each  of  the  MF  and  meteor  measurements  from  this  (unknown) 
compromise.  Thus  any  values  of  ax  and  ay  need  to  be  interpreted  both  in  terms  of 
measurement  error  and  also  the  natural  spatial  and  temporal  variability  of  the  medium 
in  which  the  measurements  are  being  undertaken,  as  well  as  the  design  of  the 
experiment  being  used  to  perform  the  measurements. 

An  alternative  way  to  consider  the  situation  is  to  consider  that  there  are  two  aspects  to 
our  measurements  -  a  component  which  is  fully  correlated  between  the  two  techniques, 
and  an  uncorrelated  component.  With  this  viewpoint,  the  fiilly  correlated  component 
represents  the  common  features  of  the  measurements,  while  the  terms  ax  and  ay 
represent  the  departures  of  the  measurements  from  this  correlated  component,  and 
therefore  are  a  measure  of  the  contribution  of  each  technique  to  the  uncorrelated 
component. 

It  is  also  apparent  that  at  some  times  of  observation  in  Figure  1  (e.g.,  the  morning  of 
June  16),  the  two  different  measurements  are  very  different  -  in  some  cases,  even 
^  opposite  in  sign.  Close  investigation  of  such  events  would  probably  reveal  some  sort  of 

large  spatial  variability,  such  as  a  large  localized  buoyancy  wave  in  the  field  of  view, 
which  might  be  seen  by  one  technique  but  not  the  other.  Of  course  such  occurrences 
must  happen,  and  deserve  further  scrunity.  However,  they  are  not  the  subject  of  our 
treatment  here,  since  we  wish  only  to  examine  statistical  differences  between  the  two 
methods  -  anomalies  like  these  points  are  simply  considered  as  part  of  the  natural 
course  of  statistical  fluctuations. 
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Despite  the  fact  that  the  techniques  being  compared  may  often  have  some  degree  of 
non-commonality,  useful  information  is  forthcoming  from  these  graphs.  First,  they 
place  limits  on  the  allowable  values  of  qq.  For  example,  if  po  =  1  is  not  an  allowable 
solution,  we  can  states  that  the  two  methods  are  different  (as  seen,  for  example,  in 
relation  to  Figure  6).  Secondly,  the  values  of  cr^  and  Uy  place  upper  limits  on  the 
intrinsic  instrumental  errors.  Thirdly,  the  relation  between  po»  and  cr^  is  clearly 
specified  by  graphs  like  Figure  5b.  The  procedure  therefore  can  be  quite  instructive  in 
relation  to  comparisons  between  different  measuring  techniques.  We  have  concentrated 
on  cases  where  the  expected  measurements  are  in  the  ratio  1 : 1 ,  or  1 1 ,  but  the 
procedure  is  quite  general.  It  can  be  used  in  many  other  cases  where  it  is  known  that  2 
measurements  are  proportional,  but  may  not  be  in  the  ratio  1 : 1 .  In  this  case,  the 
information  about  limits  on  the  optimum  values  for  the  “gain"  of  one  method  over  the 
other  (i.e.  the  value  of  po)  can  be  determined. 
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5.  The  Interpretation  of  Correlation  Coefficient 


The  correlation  coefficient,  p,  for  data  sets  {a;i}  and  {j/i}  is: 


(23) 


<  Xj  ><  pi  >  _  1 _ 

Jd+Uxi+g) 

y  X  y 


Notice  that  the  correlation  depends  on  all  of  the  parameters  ax,  ay,  go,  and  Now,  a 
quantitative  test  of  whether  two  measurements  “agree”  or  “disagree"  (e.g.,  a  statistical 
hypothesis  test)  involves  calculation  of  p  from  the  data  (using  sample  means  in  the 
numerator  and  denominator)  and  a  comparison  of  this  number  with  a  prediction  based 
on  known  values  of  ax,  ay,  go,  and  S^.  Thus  quantitative  interpretation  of  the 
numerical  value  of  p  requires  knowledge  of  all  of  these  parameters.  Without  knowledge 
of  these  parameters,  it  is  impossible  to  judge  whether  a  numerically  small  correlation 
results  because  the  instruments  are  measuring  different  physical  quantities  or  because 
the  instruments  are  measuring  the  same  quantity,  but  with  observation  errors  that 
exceed  the  variability  of  the  measured  quantity  (i.e.,  ax  >  and/or  ay  >  Sy). 

For  example,  suppose  the  observation  errors,  modelled  by  ax  and  ay,  are  equal  to  the 
geophysical  variabilities,  Eg?  and  E^,  i.e.,  suppose  ax  =  and  ax  =  E^.  In  this  case 
p  =  0.5.  Thus,  we  can  conclude  that  when  observations  errors  equal  or  exceed  the 
underlying  geophysical  variability,  measured  correlation  coefficients  will  be  smaller 
than  0.5.  This  does  not  necessarily  indicate  that  one  or  the  other  measurement  is 
“correct”  in  any  sense.  It  may  simply  indicate  that  one  or  both  of  the  data  sets  is 
“noisy”. 

For  further  illustration,  consider  the  following  “thought  problem”  wherein  the 
geophysical  variability  is  zero.  Suppose  that  the  actual  velocity  in  some  region  of  the 
mesosphere  is  exactly  zero.  Consider  two  instruments  which  make  independent 
observation  errors  and  let  those  instruments  repeatedly  probe  the  region  to  produce  a 
time  series  of  velocity  estimates.  Let  the  observation  errors  ax  and  ay  be  arbitrarily 
small  (e.g.,  suppose  they  are  on  the  order  of  1  cm/s  or  less).  Then  each  instrument  will 
produce  a  time  series  consisting  of  small  velocity  estimates,  on  the  order  of  ±1  cm/s. 
The  corrlation  between  the  two  time  series  will  be  zero  !  And  yet  each  instmment  is 
making  a  high  precision  measurement  of  the  true  velocity,  because  the  mean  value  will 
be  removed  from  the  measurements  when  computing  the  correlation  coefficient. 
Obviously,  this  is  an  extreme  example,  but  it  illustrates  that  a  small  correlation 
coefficient  is  essentially  meaningless  without  additional,  and  often  unavailable,  side 
information  about  the  geophysical  variabihty  of  the  underlying  quantity  that  is  being 
measured. 

In  summary,  high  correlation  always  indicates  agreement  between  two  techniques,  but 
low  correlation  does  not  necessarily  indicate  “poor  agreement”  between  two 
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techniques.  Thus,  correlation  coefficients  can  not  be  used  to  show  that  two  techniques 
disagree  (and  hence,  to  show  that  one  of  the  techniques  is  “wrong")  without  a  careful 
assessment  of  the  relative  magnitudes  of  observation  errors  and  geophysical  variability. 
This  often  requires  more  information  than  is  available. 
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6.  Conclusions 


A  procedure  has  been  described  which  permits  meaningful  comparison  of  different 
measurements  of  correlated  quantities,  using  a  scatter-diagram  approach.  It  is 
emphasized  that  hnear  least-squares  regression  techniques  always  result  in  an 
under-estimate  of  the  true  slope,  and  the  need  to  recognize  that  errors  exist  in  both  the 
abscissa  points  and  the  ordinate  points  must  be  recognized.  The  inter-relationship 
between  the  intrinsic  system  errors  and  the  optimum  slope  is  derived,  and  presented  in  a 
graphical  manner.  Useful  hmits  can  be  placed  on  the  errors  and  the  optimum  slope.  We 
hope  that  this  new  method  will  find  application  in  a  wide  range  of  situations,  e.g., 
ranging  from  radar  and  satellite  measurements  to  any  physical  measurements.  The 
technique  can  also  be  used  to  calibrate  multi-sensor,  radar  and  satellite,  systems. 
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